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Abstract 

We describe the computation of the amplitude A2 for a kaon to decay into two pions with 
isospin 1 = 2. The results presented in the letter fTl| from an analysis of 63 gluon configu- 
rations are updated to 146 configurations giving ReA2 = 1.381(46)stat(258)syst lO^^GeV and 
ImA2 = — 6.54(46)stat(120) systlO^^'^ GeV . ReA2 is in good agreement with the experimen- 
tal result, whereas the value of ImA2 was hitherto unknown. We are also working towards a 
direct computation of the K — )■ (;r;r)/=o amplitude Aq but, within the standard model, our re- 
sult for ImA2 can be combined with the experimental results for ReAo, ReA2 and e'/e to give 
ImAo/ReAo = —1.61(28) x 10^^ . Our result for ImA2 implies that the electroweak penguin 
(EWP) contribution to e'/e is Re(e7e)EWP = -(6.25 ±0.44stat± 1.19syst) x 10"^. 
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I. INTRODUCTION 

It was in K nn decays that both indirect f2] and direct [l3]-l6l CP-violation was first dis- 
covered and a quantitative understanding of the origin of CP- violation, both within and beyond 
the Standard Model, remains one of the principal goals of particle physics research. Lattice 
QCD provides the opportunity of computing the non-perturbative QCD effects in general and 
in hadronic CP-violating processes in particular. The evaluation of these effects in K ^ 
decays is an important element in the research programme of the RBC-UKQCD collaboration 
and in this paper we report on the evaluation of the (complex) decay amplitude A2, correspond- 
ing to the decay in which the two-pion final state has isospin 2. This is the first realistic ab 



initio calculation of a weak hadronic decay. Our final result can be found in Eq. (25 ), which we 
reproduce here for the reader's convenience: 

ReA2 = 1.381(46)stat(258)systlO-^GeV, ImAz = -6.54(46)stat(120)systl0'i3GeV. (1) 

This is an update of the result presented recently in Ref. fP] with greater statistics (146 config- 
urations compared to 63 in [IJ). More importantly, in this paper we present the details of the 
calculation and the analysis which could not be presented in the original letter [1 J. For Re A2 we 
find good agreement with the known experimental value (1.479(4) x 10^^ GeV obtained from 
decays), whereas the value of ImA2 was previously unknown. 

This is the first quantitative calculation of an amplitude for a realistic hadronic weak de- 
cay and hence extends the framework of lattice simulations into the important domain of non- 
leptonic weak decays. To reach this point has required very significant theoretical developments 
and technical progress. These are discussed in the following sections and include: 

1 . the control of nn rescattering effects and finite- volume corrections when two hadrons are 
present in the final state; 

2. the use of carefully devised boundary conditions to tune the volume so that the decay can 
be simulated at physical kinematics; 

3. the development of techniques for non-perturbative renormalization which has made it 
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possible to calculate the matrix elements of the four-quark operators in the effective 
Hamiltonian with good precision and without the use of lattice perturbation theory; 

4. the improvement of algorithms and teraflops-scale computing which has made it possible 
to perform simulations at physical quark masses. 

It has therefore required a major endeavour to control all the ingredients of the calculation to 
arrive at the final result. The systematic errors in Eq. ^ are dominated by the simple fact 
that the present calculation was performed at a single, rather large, value of the lattice spacing 
(a ~ . 14fm). With the greatly enhanced computing facilities made available to our collaboration 
and to others, the methods described in this paper can now be used at other lattice spacings to 
eliminate, or at least greatly reduce, the lattice artefacts. 

A major goal of our research programme is to calculate directly the amplitude Aq for 
K — )■ (;r;r)/=o decays, in which the final-state pions have total isospin 7 = 0, and e'/e, the 
quantity which characterises direct CP- violation in K ^ nn decays, and we reviewed the sta- 
tus of our work in [7J. The evaluation of Aq is considerably more difficult than the present 
calculation. Firstly, since the two-pion state has vacuum quantum numbers we must evaluate 
disconnected diagrams with sufficient precision. Secondly, in order to obtain physical kinemat- 
ics while avoiding the use of excited states, we must investigate alternative methods of inducing 
momentum in the final state without breaking isospin. (In the present calculation we do break 
isospin symmetry through the use of different boundary conditions on the u and d quarks, but 
circumvent the issue of mixing with 1 = states since the final state has no / = component 



because of charge conservation; this is explained in Sec. Ill ) Potential methods of improving 
the statistical precision in the calculation of disconnected-diagrams include the use of advanced 
propagator-generation techniques such as all-to-all propagators or low-mode/all-mode averag- 
ing. We are also investigating the use of G-parity boundary conditions [8J in order to achieve 
physical kinematics for decays into 7 = two-pion states. In the meantime, while we are de- 
veloping and implementing these techniques for the direct evaluation of Aq, within the standard 
model we can combine our result for ImA2 with the experimental values of RcAq, ReA2 and 
e' / £ to determine the remaining unknown quantity Im Aq, so that the values of both the complex 
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amplitudes Aq and A2 are now known (see Sec. |IIIB| ). We repeat however, that our ultimate goal 
is to compute Aq directly, and we look forward to presenting results from a realistic computation 
in the future. 

This indirect determination of Aq is also important in that it determines the 0(5%) contribu- 
tion of direct CP violation to £ SHOl. The relevance of such precision in tests of the Standard 
Model is due to the major recent improvement in the evaluation of the Bf^ parameter for which 
recent calculations have reduced the uncertainty to less than 3% []TT|. (see Sec. IIIB I. 

Since different authors use different conventions for the amplitudes we should state ours 
carefully. We define A/ (/ = 0, 2) by \/2Aj = {{nn)i \ Hw \K^) and the corresponding experimen- 
tal results are ReA2 ~ IA2I = 1.479(4) x 10"^ GeV and RcAq ~ |Ao| = 3.320(2) x 10"^ GeV. 
Expressions for the widths for — )• n^n^, Ks — ?■ n'^n^ and Ks — ?■ n'^n'^ decays in terms of 



the amplitudes are given in Eqs. ( [26] ), p4| ) and p5| ) and the surrounding discussion. 

The structure of the remainder of the paper is as follows. In the next section we present the 
details of the simulation and explain the properties of the ensembles which were used. This is 



followed in Sec. Ill by a description of our analysis together with the final results. A presentation 
of the technical details of some of the components of the analysis, including the determination 
of the systematic errors are postponed to later sections. The renormalization of the operators 



present in the effective weak Hamiltonian is described in Sec. IV and the remaining sections 
are devoted to a detailed discussion of the systematic errors. Since the matrix elements were 
calculated on a single coarse lattice, the corresponding artefacts are the largest component of 
the systematic error and we explain how we estimate them in Sec.|Vj In Sec. VI we discuss the 



errors due to partial quenching and in Sec. VII we present the remainder of the error budget. 



Finally in Sec. VIII we summarise and discuss the prospects for further work. 



II. DETAILS OF THE SIMULATION 



In this section we start with an explanation of the discrete QCD action used in our simula- 
tions (Subsec. fll A ). We then present the quark masses which we use and discuss the determi- 
nation of the lattice spacing (Subsec.[ITB|) and finally in SubseclnClwe discuss some technical 
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issues concerning the calculation of the correlation functions from which the required matrix 
elements are determined. 

A. Lattice Action 

For the quarks, we choose to use the domain wall fermion formulation [fT2l - [T4]| . This is a five 
dimensional description of QCD on a hypercubic grid, in which the fifth dimension of length Ls 
serves to separate the left- and right-handed fermion chiralities which appear as surface states 
bound to opposite four-dimensional faces of the fifth dimension. The elusive chiral symmetry 
is restored in the limit Ls—^°°. At finite the chiral symmetry is explicitly broken as the chiral 
modes can propagate across the fifth dimension. The symmetry breaking can be parametrised 
by the quantity mres, the residual mass, which additively renormalizes the bare quark masses. 
Its magnitude is governed by the density of eigenmodes of the 4D Hamiltonian obtained from 
the transfer matrix in the fifth dimension [fTSl . The contributions of the extended eigenmodes 
with eigenvalues above the mobility edge (which separates the localized low-modes from the 
extended high-modes) are dominant at small but fall exponentially as Ls is increased. In 
modem simulations with large Ls, mres is dominated by the density of near-zero eigenmodes; 
these are associated with localized and short-lived dislocations, or tears in the gauge fields 
which cause a change in the gauge field topology. 

We now discuss the choice of the gauge action. Until recently, our simulations [fT6l [TTl 
have been performed using the Iwasaki renormalization-group improved gauge action, which 
has been shown to allow adequate gauge-field topology change while retaining good chiral 
properties in Monte Carlo simulations when used in conjunction with domain wall fermions. 
The lightest (unitary) pions in these simulations had masses of about 290 MeV and the results 
were extrapolated to the physical value, m;^ ~ 140 MeV. In the present computation of ^ — nn 
decay amplitudes, we perform the simulations with sufficiently light quark masses that the pions 
have (almost) their physical masses. However as the quark masses are decreased, the pions 
propagate over larger distances and they are more strongly affected by finite-volume effects; 
this necessitates the use of physically larger lattices. In order to make the simulation affordable. 



6 



the large lattice is achieved by increasing the lattice spacing a (decreasing the inverse gauge 
coupling j8 used in the simulation). Unfortunately, as /3 is lowered, the dislocations appear 
more frequently and thus mres becomes large. To counter this effect we modify the Iwasaki 
gauge action with a weighting factor known as the Dislocation Suppressing Determinant Ratio 
(DSDR) |fT8l - l2T]| . allowing us to tune the molecular dynamics force in the gauge evolution 
to suppress configurations with large numbers of near-zero modes while retaining adequate 
topological change. This is discussed in more detail in ref. [|22|. For the remainder of this 
paper we label this action and the corresponding ensembles by IDSDR (representing Iwasaki + 
DSDR). The gauge action and ensembles without the DSDR correction are referred to simply 
by the label "Iwasaki". 



B. Parameters of the Simulation 



We have generated two ensembles of 2 + 1 flavor domain wall fermions with the IDSDR 
gauge action at /3 = 1.75 (corresponding to = 1.364GeV, see below) and a lattice size of 
32^ X 64 X 32, where the final number is L^, the length of the fifth dimension. We determine 
the residual mass to be mres = 0.001843(8), approximately equal in size to the 3.6 MeV aver- 
age of the up and down quark masses [22J. (Masses written without explicit units are to be 
understood as being in lattice units, so that for example, mres = 0.001843(8) should be read as 
^"Wres = 0.001843(8).) The ensembles are generated with a simulated strange-quark mass of 

= 0.045 and have light-quark masses of m; = 0.001 and m/ = 0.0042, with corresponding 
unitary pion masses of approximately 170 MeV and 250 MeV respectively. For the determina- 
tion of the lattice spacing a and the physical bare quark masses used in the current project, as 
well as for the computation of the particle spectrum, decay constants and the kaon bag parame- 
ter Bk, we generate quark propagators with three heavy valence masses, 0.055, 0.045 and 0.035, 
and four light valence quark masses, 0.008, 0.0042, 0.001 and 0.0001. The lightest partially- 
quenched pion has a near-physical mass of approximately 140 MeV. The analysis presented in 
this paper is performed using 146 configurations from the 0.001 ensemble, each separated by 
8 molecular dynamics time units, with additional strange quark propagators with rrih — 0.049 
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corresponding to our original estimate of the physical value of the (bare) strange quark mass, 
and light-quark propagators with a valence mass of 0.0001 . The subsequent detailed analysis 
with greater statistics and improved procedures have yielded the value 0.0472(6) for the bare 
physical strange quark mass. 

We obtain the lattice spacing and the two physical quark masses m^,^ and using a com- 
bined analysis of these IDSDR ensembles and our 32^ x 64 x 16 and 24^ x 64 x 16 domain 
wall fermion configurations with the Iwasaki gauge action at /3 = 2.25 and /3 =2.13 respec- 
tively [fT6l[T7ll . This involves a combined fit of the pion and kaon masses and decay constants 
and the mass of the f2-baryon as functions of the quark masses and lattice spacing. We use 
three different ansatze for the quark-mass dependence in order to estimate the systematic er- 
ror on the chiral extrapolations. Two of these are obtained from next-to-leading order (NLO) 
partially-quenched chiral perturbation theory with and without finite-volume corrections, and 
the third assumes a simple linear mass dependence (labelled analytic in the following). Follow- 
ing our 2010 analysis ifTTl of the two Iwasaki lattices, we extrapolate to the continuum limit 
along a family of scaling trajectories (lines of constant physics) that are defined by constant 
values of ni;;^, mx and niQ^; i.e. by imposing the condition that these masses have no lattice 
cutoff dependence on the scaling trajectory. The leading dependence on a of the remaining 
quantities is expected to be 0{a^) and in our fits we assume such a quadratic dependence. Note 
that the coefficients of the a^ terms are not constrained to be equal for the two different lat- 
tice actions. From the combined chiral and continuum fits we determine the lattice spacings 
and physical quark masses required for the pion, kaon and Q. masses to match their physical 
values, obtaining for the IDSDR ensembles an inverse-lattice spacing of = 1.364(9) GeV 
and dimensionless physical quark masses of m; = 0.00178(3) and = 0.0490(6), which cor- 
respond to 3.09 ±0.11 and 84.1 ±2.0MeV respectively when expressed in physical units in 
the MS scheme at 3 GeV. Here m = m + mres and the quoted errors contain both statistical and 
systematic contributions estimated using the procedures developed in ref. [fTTl . 

The numbers presented above were all obtained from an analysis of the 146 configurations 
used below in the evaluation of the K ^ nn matrix elements. Ref. [|22ll contains a detailed 
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analysis on an extended set of ensembles (including 180 configurations for m/ = 0.001). The 
corresponding values in Ref. Il22ll include a^^ = 1.371(8) GeV for the inverse lattice spacing, 
rhi = 0.00176(2) and = 0.0486(6) for the dimensionless physical quark masses and 3.05 ± 
0. 1 1 and 83.6 ±2.1 MeV respectively for the quark masses in physical units in the MS scheme 
at 3 GeV. 

In order to correctly propagate the correlations between the data used in the determination of 
the lattice spacings and physical quark masses with that of the present calculation of the K — > nn 
matrix elements we make use of the super-jackknife method, in which the statistical fluctuations 
associated with each ensemble are maintained separately, and the total error is determined by 
combining these contributions in quadrature. This prevents accidental correlations between the 
statistically independent data on each of the ensembles, and therefore improves on the bootstrap 
and standard jackknife methods for combining independent data. (The super-jackknife also 
does not require the number of samples on each ensemble to be the same, a limitation of the 
traditional jackknife.) A clear description of the super-jacknife technique can be found in [|23ll . 

C. Evaluation of the Correlation Functions 

We now explain some technical details concerning the evaluation of the correlation functions 
from which the matrix elements for K — )• nn decays are evaluated. Quark propagators with 
periodic and antiperiodic boundary conditions in the time direction were computed on each 
configuration with a source at f = 0. They were then combined so as to effectively double the 
time extent of the lattice. Meson correlation functions formed using the sum of the propagators 
with periodic and antiperiodic boundary conditions can be interpreted as containing forward 
propagating mesons originating at time t = 0, whereas those calculated with the difference can 
be interpreted as containing backward propagating mesons originating from a source att = 64. 
The purpose of this procedure is to suppress the so called "around the world" effects. An 
example of such effects can be seen in the two-pion correlation function, Cyzj^{t): 

Cnnit) = (0|7;,;,(04;r(0) |0) = |(0|7;,;,(0) | ;r;r)| V^-^ + ■ ■ ■ . (2) 
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FIG. 1: An illustration of around-the-world effects in the K —?■ Tin correlation function. In the left- 
hand figure the two-pion operator at t = 64 annihilates one pion and creates another, while the weak 
Hamiltonian annihilates the kaon and a pion and creates a pion. The right-hand diagram illustrates the 
K ^ Tin transition whose matrix element we evaluate. 

The term on the right-hand side of (|2]) corresponds to the creation of two pions at time zero by 
J^ji and their annihilation by J^ji at t. The corresponding functional integral however, also has 
a contribution where each of Jijt{0) and J^ni^) annihilate one pion and create another, so that 
a single pion propagates across the entire lattice. This contribution to the correlation function 
is independent of t, and although it contains the small factor e^^^^ , where T is the temporal 
size of the lattice, it may nevertheless lead to a loss of precision. Combining the propagators 
obtained with periodic and antiperiodic boundary conditions effectively replaces T by IT thus 
suppressing this unwanted contribution. A similar effect can occur in the K — )• rnz correlator if 
the weak operator in the effective Hamiltonian annihilates the kaon and one pion and creates a 
new pion, before the two-pion interpolating operator annihilates this pion and creates another 
(see Fig.[T]). Strange-quark propagators, with periodic -i- antiperiodic combinations, were gener- 
ated with sources at tx = 20, 24, 28, 32, 36, 40 and 44 in order to calculate K ^ nn correlation 
functions with kaon sources at these times, while the two-pion sources remained at either t = 
ort = 64. Thus we could achieve time separations between the kaon and two pions of 20, 24, 28 
and 32 lattice time units in two different ways which increased the statistics. These separations 
were chosen so that the signals from the kaon and two pions did not decay into noise before 
reaching the four-quark weak operator. 

We end this section with an explanation of the sources which were used for the quark prop- 
agators and hence of the operators which create and annihilate the mesons. For the propaga- 
tors of the u and s quarks, which have periodic spatial boundary conditions, we use Coulomb 
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gauge-fixed wall sources. For the J-quark on the other hand we impose antiperiodic boundary 
conditions in some spatial directions and use Coulomb gauge-fixed momentum wall sources of 
the "cosine" type 

■^p,cos (x) = COS (p^x) COS (pyy) cos (p^z) . (3) 

Here the components of momentum are given by pi = n,(;r/L) where n, is an even or odd 
integer depending on whether periodic or antiperiodic boundary conditions were imposed on 
the quark field in direction /. For our lattice, the choice ni=n2 = I and ^3 = (or permutations) 
corresponds approximately to the kinematics of a physical K — )■ nn decay. 



As explained at the beginning of Sec III we use the Wigner-Eckart theorem to relate the 
physical amplitude A2 which we wish to determine to unphysical — > n^n^ matrix elements 
which we compute directly. When studying the propagation of two mesons we use the 
same cosine source for each J-quark, which introduces cross terms in correlation functions that 
couple to two-pion states with non-zero total momentum. For illustration, consider the case 
p = (;r/L,0,0) so that the product of the sources of the two J-quarks is 

■yp,cos(xi)^p,cos(x2) = COS (jXi) COS (y-^z) 

We require the two pions to have individual momenta pi = |x and p2 = — f x (or vice-versa), 
but the first and last terms on the right hand side of Eq. Q couple to two-pion states with 
total momentum 2£ and — 2£ respectively. We eliminate the unwanted terms in the two-pion 
correlation functions by using different sinks, exp(±?;rx,/L), for the two d quarks ensuring 
that they carry equal and opposite momenta which constrains the final state to have zero total 
momentum. In the K — )■ nn correlation functions, the kaon has zero momentum and the sum 
over the spatial position of the weak operator then ensures that the two-pion final state also has 
zero total momentum. 

The advantage of using the cosine sources is that it halves the number of inversions which 
have to be performed for the J-quark. Had we used the more conventional momentum source, 

^p(x)=e'P•^ (5) 
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we would have needed to perform two separate J-quark inversions with momentum +p for 
one and — p for the other. The cosine source eliminates one of these inversions. In practice 
we only compute J-quark propagators with antiperiodic boundary conditions in or 2 spatial 
directions, corresponding to pions with ground- state momenta |p| = and |p| = \/2n/L. As 
explained above, this choice is motivated by the expectation that, with our choice of quark 
masses, |p| = y/ln/L corresponds to on-shell kinematics, i.e. that the energy of the two-pion 
state is (almost) equal to niK- 

We mention one further subtlety. As explained above, the use of antiperiodic boundary 
conditions in two spatial directions for the d quark enabled us to match the two-pion energy 
with niK- It was shown in [l24l that it is sufficient to use the antiperiodic boundary conditions 
only on the valence down anti-quarks in the ;r+ mesons, and to use periodic boundary conditions 
for the sea quarks used in the simulations. Thus we can use the gluon configurations already 
generated in which periodic boundary conditions were imposed on all the sea quarks. 



III. THE ANALYSIS 



In this section we describe the evaluation of A2. While the results presented in Eq. (25) 
towards the end of this section contain our estimates of the uncertainties, we postpone the 
detailed discussion of the determination of the systematic errors to the subsequent sections. 

The generic form of the effective Hamiltonian for K {TtK)i^2 decays is 

H,s=%V:,V,dY,CiQ]'\ (6) 
v2 i 

where Gp is the Fermi constant, Vud and Vus are CKM-matrix elements, V^^ = 0.97429, 
Vus = 0.2253 and the C, are Wilson coefficient functions. The C, contain a dependence on 
T = -yf*Vfrf/V„* y„rf = 0.0014606 - O.OOO6O4O8/, as explained below. 

The three four-quark operators contributing to the effective Hamiltonian for A/ = 3/2 decays 
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are 



2(27,1) = (^"^')^ { ("'"')^ - i.d^d^)L] + {^u\ {ujdj)L (7) 

<^(8,8) = {^d')L{{uJuJ)R-{dJdJ)R} + {s'u')L{uJdj)R (8) 
^(8,8)mix = (^"^'■)L{(i<V)/,-(jV)^} + (5V)L(MV)^ (9) 

where the superscript 3/2 denotes M = 3/2 transitions and the subscripts denote how the op- 
erators transform under the SU(3)l x SU(3)^ chiral symmetry. /, j are color labels which run 
from 1 to 3 and L,R denote left and right (e.g. {sd)t{uu)L = {sY^{\ — '}^)d) {uY^{\ — Y^)u) and 
{sd)i{uu)R = {sY^{l — 7^)<i) {uY^i^ + 7^)m) with the spinor labels contracted within each pair 
of parentheses) . 

In physical — )■ n^n^ decays the third component of isospin, 1^, changes by 1/2, AI^ =1/2. 
As proposed and first explored in [|25l l26l . it is particularly convenient to use the Wigner- 
Eckart theorem to relate the matrix elements of the operators in ([v]) - (|9]) between |^+) and 
171^71*^) States to those of the corresponding operators with AI^ = 3/2 for the unphysical process 

I + i^A/=3/2 I „+\ "\/3/ + + I „A/=3/2 I „+, 

(^+^"ieA/,=i/2l^ ) = ^(^^ ieA/,=3/2l^ )■ (10) 



On the left-hand side of Eq. ( 10 1 the operator 2^^Y/2 ^"^^ three operators in Eqs. (7 )- 
(9), whereas on the right-hand side the operators ^^li^j2 operators are \/3Qi, where i runs 
over the labels (27, 1), (8, 8) and (8, 8)niix and 

G(27,l) = {s'd')L {uJdj)L, 2(8,8) = is'd')L {uJd^)R, e(8,8)mix = {u^d')R . (11) 



y/3/2 in Eq. ( 10) is the Clebsch-Gordan factor and, neglecting violations of isospin, Eq. ( 10) is 



exact. A2 can therefore be determined by computing the matrix elements of the three operators 



in Eq. (11) and indeed it is these three matrix elements which we compute directly. For com- 



pactness of notation we suppress the labels AI and A/^ on the operators both in Eq. (11) and in 
the following. 

The use of the Wigner-Eckart theorem to replace the operators in Eqs. (|7]) - (|9]) by those in 



Eq. (11) leads to very significant practical simplifications. All the quarks participating directly 
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in A/ = 3/2 decays are valence quarks and in such cases the effects of introducing partially- 
twisted boundary conditions (for which the valence and sea quarks satisfy different boundary 
conditions) are exponentially small [24J. In particular, we assign anti-periodic boundary con- 
ditions in some directions to the valence d quarks, so that the corresponding components of 
the momenta of the final-state ;r+ mesons are {In + 1) ;r/L , where n is an integer and L is the 
spatial extent of the lattice. The volume of the lattice has been chosen so that for pions with 
momenta \/2n/L, the energy of the two-pion state E^-ji is very close to niK, the mass of the 
kaon, iriK — Ejm, corresponding to a physical decay. (The total three-momentum of the kaon 
and of the two-pion state are zero.) The most significant simplification is that the two-pion 
state is the lightest one with these boundary conditions. With periodic boundary conditions, the 
lightest two-pion state is one with each of the two pions at rest and so when computing physi- 
cal K — )■ nn amplitudes one is obliged to consider excited two-pion states fl27il . Moreover, for 
—7- Tt^Tt^ matrix elements even with anti-periodic boundary conditions on one or more of 
the quark fields, the momentum of the is Imz/L with integer n, negating the advantages de- 
scribed above. Finally we note that by using anti-periodic boundary conditions one can achieve 
the kinematics of a physical decay on a smaller lattice than with periodic boundary conditions. 

We now turn to the determination of the matrix elements. The pion and kaon two-point 
correlation functions at zero momentum are fit to the form 

Cp{t) = (0|7p(04(0) |0) = (^e--^^-f-e-'"Hr-f)j , (12) 

where T = 128 is the total effective time extent of the lattice and mp is the mass of pseudoscalar 
meson P. Jp and jj, are interpolating annihilation and creation operators for the meson P and 



Equation (12) defines Zj^ and Zk for P = n and P = K respectively. For both the pion and 
kaon the final results are obtained by fitting between t = 5 and t = 63. The masses extracted 
from these fits are superimposed on the effective mass plots in Figs. 2(a) and |2(b)[ and the 



numerical results are given in Tab.|lj The effective mass in these plots, mp eff is defined by 
C/>(0/C/>(? + 1) = cosh(mp,eff(? - 7^/2))/ cosh(mp,eff(? + 1 - 7^/2)). 

The pions in the final state for K — )■ TtTt decays have momentum |p| = y/ln/L and in Fig.js] 
we plot the effective energy for a pion with this momentum. Since the correlation functions 
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FIG. 2: Effective mass plots for the pion and kaon. Results for mj^ and niK obtained from the fits of the 



correlation functions to Eq. ( 12 ) are shown as the horizontal hnes in each plot. 



units 


nin 


rriK 






Enn2 




lattice 

MeV 


0.10421(22) 
142.11(94) 


0.37066(68) 
505.5(3.4) 


0.17386(91) 
237.1(1.8) 


0.21002(43) 
286.4(1.9) 


0.3560(23) 
485.5(4.2) 


0.0146(23) 
20.0(3.1) 



TABLE I: Results for meson masses and energies. The subscripts 0, 2 denote p = and p = \/^n/L 
respectively, where p= |p|. 



become noisier when the pion has a non-zero momentum, we now fit over the time interval 
t = [5,35] where we can neglect the contribution from the backward propagating pion and use 
the form, 

C„{t,p = y/lK/L) = \Z^{p = V27t/L)\^e-^-' , (13) 

where p = |p| and Ey^ is the corresponding energy. The value Ejj:^ = 0. 17386(91) obtained from 
the fit (see Tab.|l]) is nicely consistent with the (continuum) dispersion relation for a pion with 
mass 0.10421(22). The subscript 2 in Ej^^2 indicates that the momentum of the pion is 
i.e. that anti-periodic boundary conditions have been imposed on the d quark in two directions. 

Next we consider the two-pion correlation function which has a larger statistical error. Hav- 
ing suppressed the around-the-world contributions by combining propagators with periodic and 
antiperiodic boundary conditions in time and neglecting the contributions from excited states. 
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FIG. 3: Effective energy plot for a pion with momentum p = \/^n/L. The horizontal hne corresponds 



to value of Ej^ obtained from a fit to Eq. (13 1 



the expected behavior of the two-pion correlation function is 

C;,;,(0 = (0|7;r;r,eW4;r,c(0)|0) = ^|Z;,;,,e|' (e"^-^ + , (14) 



where the labels c and e refer to the cosine and exponential sources discussed in Sec. II C and 
Rtw is the number of directions with anti-periodic boundary conditions on the d quark. The 
leading around-the-world effects would manifest themselves as a time-independent constant on 



the right-hand side of Eq. (14). 

We find it effective in reducing the statistical errors to calculate the quotient of two-pion and 
single-pion correlators and fit the ratio to the form 



C 



(15) 



where AE = {E^^ — 2E^) and 



ytw^z V " energy difference /S.E is not equal to zero 



because of the repulsive interaction between the two pions with isospin 2 in a finite volume. 
The two-pion energy E^t^ is then given by Ej^^^ = + lE^^, and Zt^^^^ is found from 



"tw . 



Z;r;r,e — (2 2 )Z^R. 



(16) 



We can use Eq. ( [T5] ) for values of t which are sufficiently large to neglect excited states and 
sufficiently smaller than r/2 so that the backward propagating states (and the around-the-world 
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FIG. 4: The ratios C„j:{t)/{Cn{t)f defined in Eq. ([isj) at = (left-iiand plot) and at p = \/2k/L 
(riglit-liand plot). The minimum seen in the left-hand panel around t = 52 results from the different large- 
time behavior of the numerator and denominator. While the denominator decreases exponentially as t 
increases from to 64, the numerator contains a small f-independent constant (caused by one backward 
propagating pion) which lessens its decrease at large time. If examined for < ? < 128 the ratio shown 
in the left-hand panel is symmetrical about the point t = 64. 

effects) can also be neglected. In practice, in order to improve the statistical precision, we fold 
the correlation functions, averaging the equivalent results at t and T — t. We calculate the ratio 
in Eq. ([15]) for p = 0, in which case and are just the normalization factor and pion mass 



found from the fit to Eq. ( 12 1 and for p = y/ln/L in which case Zji and Ej^ are taken from the fit 



to Eq. ( 13 1. The fit regions for the quotients are f = [5,48] for p = and [5, 22] for p = y/ln/L. 
Plots of the quotients at the two values of p are shown in Figure |4j The results for all the 
meson masses and energies are presented in Tab.|l} We also present the results for mx — E^^ji to 
demonstrate that our kinematics are close to being energy conserving. 

The momentum kjz of each pion in the two-pion state is defined from the two-pion energy 
using the dispersion relation Ejm = l\Jm^-\-k^. The interactions between the two pions lead 
to kn being different from or y/ln/L. 

Next we come to the evaluation of the K ^ nn matrix elements. In the calculation as de- 
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(a) (27, 1) operator (b)(8, 8) operator (c)(8,8)mix operator 

FIG. 5: The ratios defined in Eq. ( 17 1 for p = Vln/L. The two-pion source is at f = while the kaon 



source is at = 24. The dashed line shows the error on the fit 

scribed below, we place the two-pion source at time tjtn = (or equivalently at 64) and vary 
the position of the kaon source tx- The operators of the weak Hamiltonian are inserted between 
tjijt and tK- The symmetries of lattice QCD (including translation invariance and time-reversal) 
allow us to translate the results into K — )■ nn matrix elements. 



For each of the three operators Qj in Eq. (11), where i labels the operator, the corresponding 



K —7- Tin matrix element = {n^n^ \ Qi |^+) is extracted by calculating the ratios 



(17) 



and fitting to a constant in time t. The quantity Oj(^^^ is the K — )• %% correlator with the operator 
Qi inserted at t and the kaon and two-pion interpolating operators placed at fixed times and 
respectively. and Z;j;j e are determined from the kaon and two-pion correlation functions 
using eqs. ( [T2] ) and ( [14] ). For illustration, the left-hand side of equation ( fTTj ) is plotted in Fig. [5] 
for each of the three operators for the choice = 24. The figure demonstrates that sufficiently 
far from the kaon and two-pion sources the data is indeed consistent with the expected constant 
behavior. We determine the matrix elements by fitting the data between f = 5 and ? = — 5, 
where t denotes the time distance from the two-pion source. The results for {7.^7.^%,^) 
obtained from the fits are indicated on the plot together with their errors. 

The finite-volume matrix elements computed in the lattice simulations are related to the 
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p 


Enn (MeV) 


(MeV) 


qTz 


8 (degrees) 




y/2n/L 


286.4(1.9) 
485.5(4.2) 


17.63(36) 
196.8(2.2) 


0.0659(13) 
0.7350(72) 


-0.311(18) 
-7.96(2.07) 



TABLE II: The two-pion energy E^ji, k^, qn and i'-wave phase shift 

corresponding infinite- volume ones £/i by the Lellouch-Liischer factor [|27l l28l : 

2 



I d(j) ^ d5 



2nqn V f^Qn dq,r 



(18) 



where the quantity in square brackets (denoted by LL in Tab. Ill ) contains the effects of the 
Lellouch-Liischer factor beyond the free-field normalization. 5 is the 5- wave phase shift, is 
a dimensionless quantity related to the pion momentum k^j^ by q^^ = kT^L/ln and ^ is a kinematic 
function defined in [27|. Once has been measured and qy^ determined, 5 can be calculated 
using the Liischer quantisation condition [I29M : 



nn = d{kjr)+^{q^) 
Results for E:^^, k^i, q-n and 5 are presented in Tab.|ll} 



(19) 



Since / dq-n can be calculated analytically the only unknown in equation ( fTSj ) is dS / dq-j^. 
The results for the phase shift are plotted against k^ and compared with experimental results 
[|30l[3n in the left-hand plot of Fig.[6[ we see good agreement. Near p = we assume that 5 
is linear in k^^ in order to calculate d5 / dqjt (see the right-hand plot of Fig. [6]). For p = Vln/L 
we use the phenomenological curve [[32ll shown in Fig.|6]to calculate the derivative of the phase 
shift at the corresponding value of q:^. The derivative of the phase shift is found to be a small 



term in comparison with (9^/ dqj^. Results for d(j)/ dq-n and d5 / dqj^ are presented in Tab. Ill 
We perform the analysis for four separations 5t between the kaon and two-pion sources, 

5t = 20, 24, 28 and 32. The physical decay amplitude A2 is given in terms of the matrix elements 
by 

(20) 



4' = a-'^^V,,V:,l^Qi^)Z,,i^a) -^'^ 
^ V2 ,-j 
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p 


d(j)/dqn 




LL 




y/ln/L 


0.2413(90) 
5.014(21) 


-0.0824(32) 
-0.2911(23) 


0.9632(14) 
0.9411(71) 



TABLE III: Contributions to Lellouch-Liischer factor. The second and third columns provide numerical 
values for two of the quantities entering the Lellouch-Liischer factor given within the square brackets in 



Eq. ( 18 1, while the fourth column gives the value of the complete factor. 
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FIG. 6: Plots of the 1 = 2 two-pion s-wave phase shift against momentum k^. Our results at p = and 
p = ^/2k/L are denoted by the red circles and the dashed curve is the phenomenological representa- 
tion from ref. Il32l . The left-hand plot is a comparison of the calculated phase shift with experimental 
results |[30]432l . The right-hand plot is a zoom into the small region, demonstrating the approximate 
linear behavior of the phenomenological curve in the region of /? = 0. The scattering length used in the 
straight (dotted) line is calculated using chiral perturbation theory 



where we have added the label St to indicate the separation being used and the labels / and 



i run over the three operators in Eq. (11 1. C, are the Wilson coefficients, which are generally 
calculated in schemes based on dimensional regularization; we take them to be in the MS-NDR 
scheme. The Zij are the renormalization constants which relate the bare weak operators defined 
in the lattice theory (where the lattice spacing a acts as a cut-off) to those in the MS-NDR 
scheme at scale /i. The (27, 1) operator renormalizes multiplicatively, whereas the (8,8) and 
(8, 8)niix operators mix under renormalization. The calculation of the Zij is described in detail 
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in Sec. IV and involves a non-perturbative calculation of the renormalization constants in RI- 
SMOM schemes, step-scaling to run the results to {1 = 3 GeV and matching perturbatively to the 



MS-NDR scheme at 3 GeV. As explained in Sec. IV four possible choices for the intermediate 



RI-SMOM schemes are considered. The results presented in Tab. IV are calculated using the 



renormalization constants with the intermediate scheme (ly, Iq) = ^) (see Sec. IV I. The factor 



of on the right-hand side of Eq. (20) is needed to convert from the unphysical — )• 



;r+;r+ amplitudes back to the physical — )■ n^n'^ amplitudes. 



Results for ReA2 and ImA2 for the four different separations 5t are shown in Tab. IV for 
the (almost) physical choice p = y/ln/L. Our final result for A2 is an error weighted average 
(EWA) over the four separations, defined by 



,EWA LSf^lVl^Sf 



where e§t is the statistical error in the evaluation of Af'. 



(21) 



5t 


ReAzCunits of 10"*^ GeV) 


ImA2(units of 10"" GeV) 


20 


1.411(56) 


-6.59(19) 


24 


1.346(64) 


-6.67(22) 


28 


1.427(73) 


-6.28(25) 


32 


1.295(94) 


-6.56(33) 


EWA(a) 


1.381(38) 


-6.54(15) 


EWA(b) 


1.381(44)(12) 


-6.54(19)(42) 



TABLE IV: Final results for A2. The errors on each , on EWA(a) and the first error in EWA(b) (EWA 
= error weighted average) are the statistical errors only. In the EWA(b) result the second error is that due 



from the uncertainty in the evaluation of the renormalization constants as explained in Sec. IV below. 



The errors in the results labelled by EWA(a) in Tab. IV are due to the statistical fluctuations 



on the calculated using Eq. (18). In the row marked EWA(b) the first error combines the 
uncertainty due to these fluctuations with the statistical uncertainty in the value of the lattice 
spacing and the second error is Az, which arises from the statistical uncertainty in the evaluation 
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of the renormalization constants Z/,. This is calculated using: 



^1= [^(27,1) 52(27,1) ^(27,1)]^ + ^ 



1 2 



(22) 



l2 

where /, j run over (8,8) and (8,8)niix and the dZ are the statistical uncertainties in the corre- 
sponding renormalization constants as explained in Sec.|lVj The presence of the four terms in 
the sum over / and j reflects the mixing of 2(8,8) ^i^d 2(8.8)n,ix under renormalization. ^(27,1), 



■2/(8 8) and ^(8,8)mix '^'^ right-hand side of Eq.(22) are obtained from the corresponding bare 



matrix elements using Eq. ( [T8| ). The numerical results presented here were obtained by using 
the statistical errors e^f in the evaluation of A2 so that for example: 



and similarly for the remaining operators. We have checked that performing the error weighted 
average on each operator using the statistical error corresponding to the operator makes only a 
negligible difference to the estimate of the final errors. 

For the Wilson coefficients we use the standard notation Q = z,(/i) + Tyi{iJ.) where, as ex- 
plained above, T = —VtlVtd/V*sVud- The Wilson coefficients are calculated using the equations 
in [|34ll . which uses a 10-operator basis for the effective Hamiltonian. The equations in |l34l are 
based on the pioneering Next-to-Leading Order QCD and QED calculations from the Munich 
and Rome groups [|35l - [37l . The Wilson coefficients in the 10-operator basis are related to the 
three A/ = 3/2 Wilson coefficients by 

Ci(m)+C2(m) , C9(m)+Cio(m) „ . . CtOO r (.,\ 

<~'{27,l)W = ^ + 2 ' ^(8.8)^^^ " ^(8,8)mix^A^) = ■ 

(24) 

Results for Zi and yi at /i = 3 GeV in the MS-NDR scheme are presented in Tab.fv] We ob- 
serve that the Wilson coefficients are sensitive to the value of a^. This calculation is based on 
(3 GeV) = 0.24544 which is found by solving the 4-loop running formula for [[38ll with 
initial condition ai^-* (Mz) = 0. 1 1 84 for Mz = 9 1 . 1 876 MeV [39] . The superscript (n) indicates 
the number of flavors. 
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weak operator 




Ji 


Gi 


-0.241415 





Qi 


1.11228 







-0.00392423 


0.0211096 


Qa 


0.0169695 


-0.0558734 


Qs 


-0.00349963 


0.0117843 


Qe 


0.0120747 


-0.0610235 


Qi 


0.0000940198 


-0.000161911 


Gs 


-0.000104478 


0.000652032 




0.0000275290 


-0.0103828 


Gio 


0.0000798557 


0.00243775 


2(27,1) 


0.290342 


-0.00397252 


2(8,8) 


4.70099 X 10"^ 


-8.09555x10-5 


2(8,8)„i, 


-5.22390x10"^ 


3.26016 xlO-4 



A2: 



TABLE V: Wilson coefficients at 3 GeV in the MS-NDR scheme. 
Using the procedures described above, we obtain our final results for the complex amplitude 

ReA2 = 1.381 (46)stat(258)syst 10-^ GeV, ImAj = -6.54(46)stat(120)systl0-i3GeV. (25) 

The result for ReA2 agrees well with the experimental value of 1.479(4) x 10^^ GeV obtained 
from ^+ decays and 1.573(57) x 10^^ GeV obtained from ^5 decays (the difference arises from 
the unequal u and d quark masses and from electromagnetism, two small effects not included 
in our calculation). ImA2 is unknown so that the result in Eq. ( [25] ) provides its first direct 
determination (updating the value quoted in yj). 

A detailed discussion of the determination of the systematic errors will be presented in the 
following sections. As explained in section II B[ the statistical error was obtained by analysing 
configurations each separated by 8 molecular dynamics time units. With the aim of reducing 
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ReA2 




St 


146 bins 


36 bins 


146 bins 


36 bins 


20 


1.411(56) 


1.418(52) 


-6.59(19) 


-6.55(16) 


24 


1.345(64) 


1.344(57) 


-6.67(22) 


-6.60(20) 


28 


1.427(73) 


1.411(83) 


-6.28(25) 


-6.23(29) 


32 


1.295(94) 


1.28(10) 


-6.56(33) 


-6.58(31) 


EWA(a) 


1.381(38) 


1.386(34) 


-6.54(15) 


-6.52(14) 



TABLE VI: Final results for ReA2 in units of 10"*^ GeV and ImAj in units of 10"'^ GeV. The table 



^-13 



shows a comparison between the results obtained as in Tab. IV (146 bins each with a single configu- 
rations) and those with bin-size 4 (36 bins each with 4 configurations). The errors are statistical ones 
only. 



the correlations between successive measurements, the gauge fields were shifted by 16 lat- 
tice spacings in the time direction relative to the previous configuration prior to measuring the 
quark propagators. In order to check that shifting the gauge fields is sufficient to overcome 
potential autocorrelations, we have repeated the entire analysis, including the determination of 
the physical quark masses and lattice spacings, by binning all quantities over four successive 
measurements (32 molecular dynamics time units). This is a natural choice as it matches the 
periodicity of the quark propagator measurements. The effects of the binning are completely 



negligible. For illustration we show in Tab. VI a comparison of the results for A2 obtained with 
and without the binning. 

In the remainder of the section we present the results for each of the three matrix elements 



which contribute to A2 (Sec. Ill A I and also deduce the value of the unknown quantity ImAo 
by combining our result for ImA2 with the experimental values of e'/e and other quantities 



(Sec. Ill B). In order to explain fully our conventions, we also present the explicit expressions 



for Ao, A2 and the partial widths for the K — t- TtTt decays in terms of the matrix elements. 
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A. Results for the matrix elements 

Eq. ( [25| ) contains our final results for A2 within the Standard Model. In order to facilitate 
detailed comparisons with results from future computations and to enable our results to be 
used in extensions of the Standard Model for which the Wilson coefficient functions are differ- 
ent, we now present the results for the matrix elements themselves. The results are presented 
for operators renormalized in the MS-NDR scheme at a renormalization scale of 3 GeV. Our 
convention is that V2A2 = (^{n7t)jZ}Q\Hw\K^'^ . With this definition IA2I = -sJ^IAj^qI, where 
A+o = {7t'^7r^\Hw\K^) and the corresponding partial width is given by 

nK+^n+n^) = ^\A^o\'^. (26) 



ml. 



where 




P^'=\^- 2 + Ami, ' 



1. —7- 71^71+ matrix elements 



We start with the results for the — )■ ;r+;r+ matrix elements of the operators defined in 



Eq. ( 1 1 ) in terms of which A2 is given by 



A2 = £CK3 GeV)^^s-NDRp G^v) (28) 

where j^.ms-ndR = ( | q. and the label i runs over (27,1), (8,8) and (8,8)mix . The .<^i 
take the values 

1)^^(3 GeV) = 0.03071(97) GeV^ (29a) 
i^(f|)^°'^(3 GeV) = 0.583(33) GeV^ (29b) 
"^mT^^ G^V) = 2-64(15) GeV3 . (29c) 
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2. — )• 71+ tt" matrix elements 



Alternatively we may express A2 in terms of the matrix elements for the physical 
;r+;r° decay. In this case 

Ai = ^K,C-^£C,(3 GeV)^'p-^°^(3 GeV). (30) 

where the two-pion final state is symmetrised {^{{k'^ {p)K^ {—p) \ + {K^{—p)n^{p)\). We 
find the matrix elements to be 

■s^'^7."i)°^(3 GeV) = 0.0461(14) GeV3 (31a) 
^,MS-NDRp GeV) = 0.874(49) GeV^ (31b) 
-^Sf (3 GeV) = 3.96(23) GeV^ . (31c) 



3. Contributions to A2 from the Matrix Elements 

Finally we present the separate contributions to A2 in Eq. ( [25] ) from the matrix elements of 
the three different operators: 

Re(A2)(27,i) = (1.398 ±0.044) 10-^ GeV; Im(A2)(27,i) = (1.55 ±0.36) lO^^^ QeV 

Re(A2)(8,8) = (4.29 ±0.24) 10-" GeV; Im(A2)(8,8) = (4.47 ±0.25) lO^i^GeV 

Re(A2)(8,8)_ = (-2.14±0.12)10-iOGeV; Im(A2)(8.8)_ = (-8.14±0.47) lO-^^GeV. 

(32) 

B. Prediction for ImAo 

Before describing our indirect determination of the unknown quantity ImAg, we summarise 
our conventions. Aq is defined by a/2Ao = (^{miyi^^Q\Hy/\K^^ . Defining the amplitudes A_|_„ 
and Aqo by 

A+^ = {n+n-\Hw\Ks) and Aqq = {k^k^\Hw\Ks) , (33) 
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the corresponding partial widths are given by 



Y{Ks^TZ^n-) = ^ |A+_|2 (34) 



1 

1671 



r(i^5->^V) = T^lAool'^, (35) 



where the relative momenta are given by 



and Aqq are given in terms of Aq and A2 by 



A+_ = yl^A2e'^ + -j=Aoe'^ (37) 
Aoo = 2^A2e-^-^Aoe'^\ (38) 



where 5j is the s-wave nn phase shift for isospin /. With these definitions we now evaluate 
ImAo. 

Having obtained A2, and recalling that RcAq is known from experiment, the remaining un- 
known quantity is ImAo. We now determine this by combining our result for ImA2/ReA2 from 



Tab. VII , with the experimental values of 



Re 



£'\ CO rimA2 ImAo 



^/2\£\ [ReA2 ReAo 
(D, |e| and RcAq, where co = ReA2/ReAo and 



(39) 



where 77,- y = — — -— . (40) 

3 A{Ks n'TtJ) 



The numerical values which we use for these quantities are given in Tab. VII The systematic 
error on ImA2/ReA2 is found by combining in quadrature the systematic error on ReA2 and 
ImA2 with the error due to lattice artefacts excluded. We then add in quadrature a single estimate 
of 5% systematic error on ImA2/ReA2 due to lattice artefacts. This estimate is based on the 
Symanzik theory of improvement which implies that the artefacts are proportional to and in 
the absence of any knowledge of the constant of proportionality we use the spread of the derived 
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values of the lattice spacing in Tab VIII below as a guide. Our result and error for ImAo/ReAo 
are very insensitive to the estimate of the artefacts in ImA2/ReA2. 



Rearranging Eq. p9) we determine the unknown quantity ImAo within the Standard Model, 
finding 

ImAo = -5.34(62)stat(68)syst x 10'^^ GeV. (41) 



The error on ImAg is obtained by combining the errors on the quantities in Tab. VII in quadra- 
ture. The relative contribution to ImAo/ReAo from ImA2/ReA2 and the term containing the 
experimentally known contributions is given by: 

ImAo ImA2 Vllele' 



ReAo ReA2 (o 



-1.61(19)stat(20)systx 10-4 = -4.42(31)stat(89)syst X 10-5 - 1.16(18) X 10- 



(42) 



Thus we see that while the error on the determination of ImAo is dominated by the uncertainty 
in the experimental value of e'/e, the contribution of ImA2/ReA2 is significant (about 25%). 
Of course our ultimate aim is to calculate Ao directly and we hope to be able to report on this 
soon; an important step towards this goal was presented in 

The ratio ImAo/ReAo allows us to determine the effect of direct CP- violation in — )■ nn on 
£, customarily denoted by fCg dH, (K'e)abs = 0.924 ±0.006. where the subscript "abs" denotes 
that at present only the absorptive long-distance contribution (Im ri2) is included ifTOll (the 
error is now dominated by the experimental uncertainty in e' /£). The analogous contribution 
from the dispersive part (Im M12) [fTOll is yet to be determined in lattice QCD, but we describe 
progress towards being able to do this in [|40ll . 



Using our value of ImA2 in Eq. (25) and taking the experimental value given above 
for ReA2 from decays we obtain the electroweak penguin (EWP) contribution to e'/e, 
Re(e7e)EWP = —(6.25 ±0.44stat± l-19syst) x 10^"^ (the experimental value for the complete 
Re(e7£) is 1.65(26) x 10-^ JSl). Even though we have labelled this contribution EWP, 
and indeed it is dominated by the matrix element of the EWP operator 2(8,8)n,ix' result 



contains contributions from all three components to ImA2 in Eq. (32). The (renormalization- 
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e'/e 


(1.65 ±0.26) X 10-3 


(0 


0.04454(12) 


|£| 


(2.228 ±0.011) X 10^3 


ReAo 


3.3201(18) X 10-' GeV 


ImAi/ReAa (lattice) 


-4.42(31)stat(89)systXl0-5 



TABLE VII: Experimental values of the quantities in Eq. ([39]l which is used in the determination of 
ImAo, together with the result for ImA2/ReA2 from this paper. 

group invariant) sum of the contributions from the two EWP operators 2(8.8) and 2(8, s)^;^ is 
-(7.34±0.52stat±1.39syst) x lO^^. 

We end this section with a brief comparison of an earlier result obtained using finite-energy 
sum rules dHJ, where the contribution to e'/e from the operator 2(8,8)n,ix (renormalized at 
2 GeV) was found to be —(11.0 ±3.6) x 10^^. Our result for this particular contribution is 
(—7.88 ± 0.43) X 10^^. (Note that the contribution from 2(8.8),nix itself is not renormaliza- 
tion group invariant.) We also note that our result is consistent with expectations based on the 
vacuum saturation approximation at scales around 2 GeV [|4ni42l . For a comprehensive general 
recent review on kaon decays we refer the reader to Ref. [|43ll . 



IV. RENORMALIZATION OF THE LATTICE OPERATORS 



We have seen in Sec. Ill above, that in order to determine the physical amplitudes we need 
to combine the K — t- tztz matrix elements with Wilson coefficient functions. The coefficient 
functions are calculated in perturbation theory and most often correspond to renormalization 
schemes based on dimensional regularization, such as the MS-NDR scheme. We therefore 
need to determine the matrix elements of the weak operators also renormalized in the MS-NDR 
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scheme and schematically this is done as follows: 



bare 




operators renormalized in 




operators 


< lattice 




intermediate scheme (s) 


^ Pert.Th. ^ 


renormalized in > 


operators 




(RI-MOM, RI-SMOM) 




MS - NDR scheme. 



In the first step we perform non-perturbative renormalization (NPR) on the bare lattice operators 
to obtain operators defined in a renormalization scheme which can be simulated numerically, 
such as the RI-MOM or RI-SMOM schemes discussed below p4] - l46ll . Since we cannot perform 
simulations in a non-integer number of space-time dimensions, the introduction of intermediate 
schemes is necessary. In the second step continuum perturbation theory is used to relate the 
operators in these intermediate schemes to the MS-NDR scheme. In this way we avoid the 
use of lattice perturbation theory, which frequently converges more slowly than its continuum 
counterpart and for which it is more difficult to calculate the higher-order corrections. 

Of course the relations between the bare lattice operators and those renormalized in the 
MS-NDR scheme are, in principle, independent of the choice of the intermediate scheme. In 
practice, in addition to the remaining lattice systematic uncertainties, the fact that the matching 
between the operators in the intermediate schemes and MS-NDR is performed only at a rela- 
tively low order of perturbation theory means that there is a small dependence on the choice of 
intermediate scheme. As explained in the following subsections, we find it useful to use a num- 
ber of intermediate schemes and to use the spread of results as an indication of the uncertainties 
due to the truncation of perturbation theory. 

A. The Intermediate Renormalization Schemes 

The intermediate renormalization schemes we use are natural extensions of those we intro- 
duced in our recent study of the parameter of neutral kaon mixing Wfl . These in turn were 
based on the schemes we had introduced for quark bilinear operators in which there are no 
exceptional channels, i.e. no channels with small or zero momenta [|45l |46ll . By imposing the 
renormalization conditions on quark and gluon Green functions with no exceptional channels 
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z, a 




FIG. 7: Schematic diagram illustrating the process in Eq. (43 1. In the diagram the arrows refer both 



to the flow of the indicated flavor quantum number and also to the indicated momentum. The spinor 
(a, j8 , 7, 5) and color (/, j, k, I) labels are also indicated. 



we suppress the systematic errors due to the breaking of chiral symmetry by infrared effects. 
We now explicitly explain the schemes which we use. For all the operators we introduce two 
ways of treating the vertex renormalization and two ways of defining the wave function renor- 
malization, leading to four renormalization schemes for the operators themselves. 



The three operators which we need to renormalize are defined in Eq. (11). 2(27, i) renormal- 
izes multiplicatively, whereas the two electroweak penguin operators 2(8,8) ^i^d 2(8.8)mix rni^ 
under renormalization so that there is a corresponding 2x2 matrix of renormalization constants. 
We start with a discussion of the renormalization of 2(27, i) for which we compute the Green 
function for the process 

d{pi)s{-p2) ^ d{-pi)u{p2) (43) 

with p\ = P2 = (pi —PiY = for ^ variety of momenta satisfying this condition. The process 
is illustrated in the diagram of Fig|7]and /i is taken to be the renormalization scale. 

Let A^^'^^''^'*^'(/7i,/)2) be the amputated Landau-gauge Green function of the bare lattice 
2(27,1)' where a,/3, 7 and 5 are the spinor labels corresponding to the incoming d and 5 quarks 
and outgoing d and u quarks respectively and /, j, k, I are the corresponding color labels. Since 
2(27,1) is multiplicatively renormalizable, the relation between the bare lattice and renormalized 
operator is of the form: 

e(lvilq) y(Iv,Iq) ^(latt) ( AA\ 

(27,1) ~ ^(27,1) ^(27,1)' ^^^> 
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where Iv labels the choice of the intermediate (one-particle irreducible) vertex renormalization 
scheme and Iq the intermediate scheme for the wave function renormalization. The index "latt" 
reminds us that the operator on the right-hand side is the bare lattice operator. The overall 
renormalization constant is obtained by evaluating a trace of A with a projection operator P^^"^ 

7(Iv,Iq) _ 7(Iq)2 1 (A'^'\ 

^(27,1) {h)ij,kl (27,1) ij,kr ^^^^ 

where Z^'^^ is the wave function renormalization constant which will be discussed below. The 
two choices we make for the projection operators are labelled by ly = or ly =^ P31 : 

\f)pair)sr+ (rV)/3a(rV)5r] (46) 
'U)pa{^)5r+Uy')paU^)5r\ , (47) 



pir)ij,kl ^ 1 

128A^(A^+1; 

{li)ij,kl _ I 



where A'^ = 3 is the number of colors. These projectors are constructed to give 1 when contracted 

with the tree-level results for A^fl'il''''^'. 

ap,Yo 

For the wave function renormalization we use the schemes defined as RI-SMOM and RI- 
SMOMy^ in ref. [1461 . which for compactness of notation, we label as Iq =</ and Iq = 7^ respec- 
tively. The corresponding renormalization constants are defined as 

4^^ = ^Tr[A(;^ and Z^'^^ = j^Tr[A^yf] , (48) 

where Ay is the amputated Green function of the conserved vector current. This completes the 
description of the determination of the renormalization constant for 2(27, 1) in the four schemes 
in which each of Iq and ly are either ^ or 7^. 

We now turn to the determination of the renormalization constants of the electroweak pen- 
guin operators Qj = 2(8,8) and = 2(8,8)mix' where the notation Qj and 28 is another standard 
one and will prove convenient in the following discussion. In this case we define two projection 
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operators Pj^""^ and /g^"'* for each scheme (ly = Jfi or 



^7 


ij,kl 
aj3,y5 


piry 

^8 


ij,kl 








ij,k! 


'puy 






ij,kl 


p(#)" 









(49) 
(50) 
(51) 
(52) 

For each scheme, let Mah {a, b = 1,S) be the matrix obtained by tracing with the amputated 
Green function A" over spinor and color indices: 



Mab = [Ph]ap[ys ^^"^ap,rS 



(53) 



with an implicit sum over all repeated indices and we have suppressed the index ly = or ^ 
defining the renormalization scheme. The matrix of renormalization constants Zab {a,b = 7, 8) 
is defined by 



1 



ZM^M 







(54) 



where the matrix is the free-field expression for M. 

With a single choice of boundary conditions, the components of momenta are quantized in 
steps of In/L, where L is the spatial extent of the lattice. In order to study the momentum 
dependence of the Green functions from which the renormalization constants are calculated 
we need to take a range of values for each component of momentum. The presence of lattice 
artefacts which are not invariant under the 0(4) group (but which are invariant under the lat- 
tice discrete symmetry group) leads to irregularities in the computed momentum dependence. 
Examples of such contributions are terms proportional to pfi)/ (E/x pji)- Such terms are 

not proportional to cP'p^ (where = ll^pji) introduce a scatter in Green functions when 
plotted as functions of p^, making it difficult to extrapolate the results to the continuum limit. 
The use of partially twisted boundary conditions [l24l l48ll allows us to scale the components of 
the momenta (almost) continuously, so that (L^/'^)/(E/iP^) ^^d p^ scale in the same way and 
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the scatter is eliminated. This technique was used in our recent calculation of the Bk parame- 
ter [|47ll where it is described in detail and it is used throughout the present calculation of the 
renormalization constants. 



B. Step Scaling 

In the preceding subsection we described how we obtain the renormalized operators 
2(2/1) (j")' Qls'sj'V) and of^'^f . (ju) on the coarse IDSDR lattice, where the renormaliza- 



tion scale IJL = p\ = P2 = {p\ — Pi) (see the discussion around Eq. (43 1). In order to hmit the 
lattice artefacts on this coarse lattice (a ~ 0.14fm) /i should not be very large. On the other 
hand if we choose /i to be too small then perturbation theory cannot be used reliably to re- 
late the operators in the intermediate schemes to those in the MS-NDR scheme. The use of step 
scaling [|49ll50ll . and in particular its recent generalization to the RI-SMOM schemes being used 
in this work [|47l [STl [52ll . overcomes this last limitation as explained below. This step scaling 
approach can also be generalised to operators which mix under renormalization [|53] [54ll and 
this is applied in our calculation. 



Imagine that we use the procedure of Subsec. IV A to obtain the renormalization constants 



Z^^^^^'^^{11q) and Z^^^'^°'\iiq) , {a,b = 7,8), on the IDSDR lattice for a renormalization scale /io 
which is sufficiently small that lattice artefacts can be neglected and which is therefore likely 
to be outside of the perturbative regime. We then repeat the same renormalization procedure 
to obtain the corresponding renormalization constants, and hence the corresponding operators, 
on the finer Iwasaki lattices mentioned in Sec. [11} (Renormalization constants on the Iwasaki 
ensembles were presented in [l53l.) The benefit of doing this is that on the finer lattices we 
can run the renormalization constants non-perturbatively from /io to a larger scale /i at which 
perturbation theory can be applied. Taking Q(2i.\) as an example, we define a step scaling 
function on the finer lattices: 



^(2V J (f^^ ^' (^(^Tj)* (A^o, a, m 



(55) 

where m is the quark mass. Since we have results at two different lattice spacings on the finer 
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Iwasaki lattices we can perform the continuum extrapolation and define the continuum step 
scaling functions as 

(T(^2';;^)^(M, Mo) = limE|^;''j^j(M,Mo,«) ■ (56) 

The step scaling function 0'(27,i) (M? Mo) describes the continuum non-perturbative running of the 
4 quark operator 2(27, i) in ^ given scheme. Because it does not depend on the lattice action, we 
can use it to run the Z factor obtained from the IDSDR lattice at a low scale jio to a higher energy 
jU where perturbation theory is more convergent. Finally, the operator 2(27, i) renormalized in 
the intermediate scheme (ly, Iq) at a perturbative scale /i is related to the IDSDR lattice operator 
by: 

2(27^^) = ag-'r)V,Mo)z;5-\^)Vo) 2SS;l)- (57) 
Having obtained the operator renormalized in the intermediate schemes at a perturbative 
renormalization scale, we convert it to the MS-NDR scheme using one-loop perturbation theory 

2gf.i)(M) = ^vJ'^^lM) . (58) 

The expressions for the conversion factors ^[^-^^{^'''(M) can be found in ref. [l55l . Since these 
are known to 0{as) the determinations of 2^7 i)(m) via different intermediate schemes (Iv,Iq) 
will differ from one another at 0{aj). The difference of results calculated via different inter- 
mediate schemes provides an estimate for the size of this effect. 

For the electroweak operators the above equations become 2x2 matrix equations with the 
constants zf^^'^^'' replaced by the matrices z*^^'^'''' and similarly for the step scaling factors. 



C. Numerical Evaluation of the Renormalization Constants 



We now present the numerical results for the conversion matrices that relate our bare lattice 
operators, 2f^"'' (' = (27, 1), (8, 8), (8, 8)niix)5 to those renormalized in the MS-NDR scheme at 
the renormalization scale ji, 2f ^ (m), 

2^(m) = 5(^-i^)^^(Ai) a(^-i^V,Mo) z(i-i^)(Mo) 2^""^ (59) 
^zfff{ii)Q^'^''K (60) 
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As explained in Sec. 



IV B 



the conversion matrix [Z^^^^^^^^^ {l^)]ab is a product of the three factors 
explicitly exhibited in Eq. ( 59 1. We have studied the four different intermediate schemes (ly, Iq) 



introduced in Sec. IV A in order to estimate the uncertainty from perturbative truncation errors 
in the (Iv,Iq) to MS matching factors S^^^'^'i^^^^ (jl) , which are known at one loop [|55ll . and 
also the uncertainty from discretisation effects. 

In the evaluation of WJ\ . our study of the renormalization of the (27,1) operator 
concluded that of the four choices of (Iv,Iq) it was the non-perturbative running functions 
a^2^'^')(MiMo) and a^^^'J^ {jl.jlo) which were best approximated by perturbation theory for 
jU ~ 3 GeV. These two intermediate schemes were then chosen for the determination of the 
matrix elements of and in estimating the truncation uncertainty. In the current work, we 

again find that the running functions in the (^, and (7^, y^) schemes, now 3x3 matrices, are 
generally well described by perturbation theory. We therefore choose to adopt the same proce- 
dure as in [|47l : we take the results from the (^,^) intermediate scheme as our central values for 
ReA2 and ImA2, and use the difference between these and the results obtained in the {y^ ^Y^) 
scheme as an estimate of the uncertainty. 

In order to minimise discretisation effects in the calculation of the Z-factors on the IDSDR 
lattices, where a is large and only one lattice spacing is available, we take as a matching point 
the low scale jUq = 1.136 GeV. We obtain: 



0.443(1) 





(61) 




0.505(1) -0.114(1) 
-0.022(3) 0.231 (2)y 

^ 

0.510(2) -0.116(1) 
-0.077(6) 0.305(4)^ 

where the quoted errors are statistical only. Here and in the remainder of this section we estimate 
and propagate the statical errors by using 100 bootstrap samples. 



Vo) 



0.489(1) 





(62) 



The block structure of the matrices given in Eqs.(61) and (62) is justified by the short- 
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distance chiral symmetry of the DWF formulation which impUes that changes in action and 
lattice spacing can be compensated by multiplicative renormalization of 2 (27,1) and mixing 
between Q-j and gg. The method described above determines the five elements of the 3 x 3 
matrix Z which are expected to be non-zero and the remaining four are set to zero because the 
chiral symmetry of the theory implies that operators with different chiraUty do not mix under 
renormalization. 

The renormaUzation constants at jUq are converted to the higher scale jU = 3 GeV using step- 
scaling functions calculated on the Iwasaki lattices, extrapolated to the continuum Umit. When 
performing the continuum extrapolation, we match scales on the different lattices by interpo- 
lating the simulated data, which are very smooth on account of our use of twisted boundary 
conditions. Twisted boundary conditions also ensure that the data lie along a continuum tra- 
jectory, and with two Iwasaki ensembles we can attempt to remove 0{a^) artefacts using a 
straight-line fit. Since we have only two lattice spacings, we choose to quote a conservative sys- 
tematic error: the difference between the results on our finest Iwasaki lattice and those obtained 
by extrapolating to the continuum. In this way we obtain 



(7(^'^)(3GeV,jUo) = 



^.942 (4) (1) ^ 

0.964 (9) (13) 0.386 (20) (79) (63) 
y 0.038 (23)(17) 2.210(76)(103)^ 

^0.876 (7) (9) ^ 

0.973 (11)(6) 0.309 (16)(67) • (64) 
y 0.166(38)(50) 1.884 (84) (45)^ 

The first quoted errors are statistical, while the second are the systematic ones from the contin- 
uum extrapolation. 

The matching to the MS-NDR scheme is performed at the scale ji = 3 GeV where per- 
turbation theory is more convergent than at the conventional scale of jU = 2 GeV. Using 



(y(^,^)(3GeV,/io) 
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a^^{3 GeV) = 0.24544, we obtain for the matching factors: 

/ 

sir, ^3 GcY) = 



1.00414 
1.00084 -0.00253 
-0.03152 1.08781 



(65) 



>MS 



(3 GeV) 



0.99112 
1.00084 -0.00253 
-0.01199 1.02921 



(66) 



Multiplying these results together and propagating the systematic errors in quadrature gives our 
final result: 



^MS,(latt) 

(3GeV) 



'^0.419(2)(1) ^ 

0.479 (5) (8) -0.022 (5) (20) (67) 

y -0.047(13)(11) 0.552(19)(28)^ 

/o.424(4)(4) ^ 

0.472 (6)(8) -0.020 (5)(21) . (68) 

\^ -0.067 (23) (30) 0.572 (28) (20)^ 

For each result, the first quoted error is statistical errors, while the second is the systematic 



Z^,;^"""(3GeV) 



uncertainty due to the continuum extrapolation in (63 1. 



D. Is 3 GeV a sufficiently large momentum for perturbative matching? 

In the previous subsections we described how we calculate the renormalization constants 
relating the bare lattice operators on the IDSDR lattices to those renormalized in the RI-SMOM 
schemes at a renormalization scale of 3 GeV. This calculation is entirely non-perturbative. In 
order to obtain the physical amplitude A2 the matrix elements of these renormalized operators 
have to be combined with the Wilson coefficient functions which are calculated in perturbation 
theory, most often in schemes based on dimensional regularisation. We therefore convert our 
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results to the MS-NDR scheme at /i = 3 GeV and, since we cannot perfomi simulations in 
a non-integer number of dimensions, this conversion has necessarily to be performed using 
(continuum) perturbation theory. At present we know the conversion factor to one-loop order 



and the difference of the results in Eqs.(67) and (68) provides an estimate of the systematic 
error due to the truncation of the perturbative matching to one-loop order in going from the 
RI-SMOM to the MS-NDR schemes. In this subsection we investigate further whether 3 GeV 
is a sufficiently large scale at which to use perturbation theory. We do this in two ways. Firstly 
we study how well the non-perturbative running tracks perturbation theory in the vicinity of 
jU = 3 GeV. We then check whether the infra-red chiral symmetry breaking effects are small at 
3 GeV. 



1. Comparing perturbative and non-pertubative running. 

It is instructive to start with the four plots in Fig. [8] These represent the running of the step 
scaling functions for 2(2?, i), determined non-perturbatively, normalized by the LO or NLO 
perturbative expressions for the four RI-SMOM schemes considered in this paper. The ratios 
are fixed to be 1 at /i = 3 GeV where we match perturbatively to the MS-NDR scheme. We 
see that for the ( ^, ^) scheme the running is very much as expected from NLO perturbation 
theory (and indeed LO perturbation theory) in the vicinity of 3 GeV and this was the primary 
reason why our central values for Bk (which is also obtained from the matrix element of an 
operator which transforms as an SU(3)lxSU(3)^ (27,1) and is related to the operator studied 
here by a chiral rotation) were quoted using ( ^, ^) as the intermediate scheme [|47l . The (7, 7) 
scheme shows a reasonable agreement between the perturbative and non-perturbative running 
in the vicinity of 3 GeV and we used the results with this intermediate scheme to estimate the 
truncation error of the matching to the MS-NDR scheme. 

For the electroweak penguin operators, while the numerical details are different, the same 
general features are also present. For illustration we present the results for the diagonal terms, 
which are the most important ones, in Figs.|9] and 10 and we follow the same procedure in 



quoting our central values and systematic errors. More details can be found in ref. [l54ll . 
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FIG. 8: The running of the step-scaUng function divided by the LO or NLO perturbative expression 
for the operator Q(2i,\) for the four RI-SMOM schemes considered in this paper. The ratio is set to 1 at 
jJ. = 3 GeV. The non-perturbative results have been extrapolated to both the chiral and continuum limits. 

2. Infrared chiral symmetry breaking effects 



We initially impose the RI-SMOM renormalization conditions at the relatively low scale of 
/Xo = 1-136 GeV where we might expect that infrared effects due to the spontaneous breaking of 
chiral symmetry may not be negligible, even after the quark masses are set to zero. This does not 
matter however, since we do not need to introduce perturbation theory until we have run all our 
results to 3 GeV. Recall that we have determined the renormalization constants needed to relate 
the bare lattice operators defined on the coarse IDSDR lattices to the RI-SMOM renormalized 
operators completely non-perturbatively, including the infrared effects. At 3 GeV we would 
expect that these effects are very small, indeed we require this to be the case in order later to 
apply perturbation theory at this scale. To illustrate that this is indeed the case, we study the 
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FIG. 9: The running of the step-scaling function divided by the LO or NLO perturbative expression for 
the 77 element of the step scaling function for the four RI-SMOM schemes considered in this paper. The 
ratio is set to 1 at = 3 GeV. The non-perturbative results have been extrapolated to both the chiral and 
continuum limits. 

size of the "wrong chirality traces" as we now describe. 

For the purposes of this discussion it is convenient to modify the projection operators defined 
in Sec. |IV A so that the tree-level projections give the identity. The 5 renormalization conditions 



imposed in Sec. IV A can be written in the schematic form: 



^(27,i)Af27,i)(M)=^ and PjAf{ii)=G,j where /, 7 = 7,8. 



(69) 



The A^ are the Green functions of the operators renormalized in one of the RI-SMOM schemes 
at the renormalization scale jU (the superscript R stands for Renormalized), the P are the projec- 



tors as defined in Sec. IV A and the constant F and constant 2x2 matrix G correspond to the 



tree-level values of the traces (with the normalization factor in Eqs. (46) and (47) the constant 
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FIG. 10: The running of the step-scaUng function divided by the LO or NLO perturbative expression 
for the 88 element of the step scaling function for the four Rl-SMOM schemes considered in this paper. 
The ratio is set to 1 at /i = 3 GeV. The non-perturbative results have been extrapolated to both the chiral 
and continuum limits. 

F = 1, but we leave the value unspecified for this general discussion). We now modify the 
projectors to 



1 



^27,1) = ;^^ ^(27,1) and P] = {PG-')j for J =7,8, 



(70) 



in terms of which the conditions ( 69 1 on the Green functions read 



^(27,i)^(27,i)(A^) = 1 and P]Af{^) = 5ij v^herciJ = 7,S. (71) 
We now introduce the 3 x 3 matrix A/,j(/i), where the labels /, j = 1,2,3 correspond to the three 
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operators (27,1), 7 and 8 respectively; 



A 



1 fl(Ai) b{^) 
c{il) 1 . (72) 
\d{fi) 1 ^ 

At 3GeV, as explained above, we require the wrong chirality constants a{3GeM), b{3GeV), 
c(3GeV) and (i(3GeV) to be small and indeed this is what we find. For example, in our pre- 
ferred ( ^, ^) scheme on the 32^ Iwasaki lattice in the chiral limit we obtain 

^ 1 -2 (2) X 10-5 2 (2) X 10-5^ 

0(2)xl0-5 1 

y-4(2)xl0-5 1 J 

Had the wrong-chiraUty traces not been small, we would have expected similar infrared effects 
in the renormalization conditions themselves and not been able to apply perturbation theory at 

this scale. 

At the lower scale of /io = 1.136GeV we expect the wrong chirality traces to be larger and 
this is indeed the case, although we find that they are actually still small. The key point here is 
that they are physical and therefore should be the same for all lattices. We find 



M323 (3GeV) 



(73) 



MiDSDR(l-136GeV) 



1 -0.002(2) -0.004(2) 
0.002 (2) 1 
-0.002(4) 1 



(74) 



for the IDSDR lattices and 



M243(1.136GeV) = 



M323(1.136GeV) 



1 -0.001(1) -0.007(1) 
0.001 (1) 1 
\^ -0.004 (2) 1 

^ 1 0.000(1) -0.006(2)^ 

0.003(1) 1 
\^ -0.008 (3) 1 J 



and 



(75) 



(76) 
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Quantity 


ChPTFV 


Analytic 


ma 


1.364(8) GeV 


1.362(11) GeV 


fn 


1.410(27) GeV 


1.386(19) GeV 


fK 


1.413(29) GeV 


1.392(28) GeV 


n) 


1.357(4) GeV 


1.362(7) GeV 



TABLE VIII: Values of the inverse lattice spacing obtained using different physical quantities 
to set the scale. For the Sommer scale ro we use the value ro = 2.433(50)(18)(13)GeV^^ = 
0.4795(99) (35) (26) fm from our detailed analysis in ll22l . The two columns of results correspond to 
the use of finite-volume SU(2) chiral perturbation theory and the analytic ansatz for the light-quark mass 
dependence. 

for the two Iwasaki lattices. Within the errors, the results are indeed consistent with our expec- 
tations. 



V. ESTIMATING THE ERROR DUE TO LATTICE ARTEFACTS 

We now begin a detailed examination of the systematic uncertainties leading to the estimates 



in Tab. XI In this section we study the largest single contribution to the systematic uncertainty, 
that due to the artefacts. 

Our calculations of the K — t- tztz amplitudes were performed at a single, rather large, value 
of the lattice spacing, a^^ = 1.364(9) GeV. As described earlier, this value of the lattice spacing 
was obtained in our standard way using the mass of the Q-baryon to set the scale and the 
masses of the pion and kaon to determine the physical quark masses. With the action which we 
are using, all other computed physical quantities have errors of 0{a^), but without a simulation 
at a second lattice spacing we cannot determine these lattice artefacts directly. In this section 
we describe our indirect estimates of the 0{a^) effects. 

We use two (related) methods to estimate the artefacts. In the first of these we imagine using 
quantities other than m^j to set the scale and observe the corresponding variation which we 
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ascribe to artefacts. The results are presented in Tab. VIII The difference between the largest 
and smallest entry in the table is about 4%. Recalling that the K — )■ tztz matrix elements are of 
dimension 3, we would estimate the corresponding uncertainty in the amplitudes to be 10-15%. 
On the other hand, it could be argued that we don't know the physical value of tq very well and 
that we should simply impose that we obtain the same value of tq on the Iwasaki and IDSDR 
lattices. This then fixes the ratio of lattice spacings on the two ensembles. Combining this 
ratio with the well determined lattice spacing on the Iwasaki ensembles from mn leads to the 
IDSDR value = 1.363(22) GeV, closer to those obtained from ma and the decay constants. 
Although this may suggest that the 10-15% estimate is conservative, because of the indirect 
nature of these estimates, we prefer to be conservative when quoting the uncertainties. 

As a second approach we set the scale from as usual and study the matrix element 
^A5=2 ^ I^K\sy^{l - f)d) {sf'il - f)d)\K^) on the Iwasaki and IDSDR lattices. This ma- 
trix element gives the dominant contribution to the indirect CP-violation parameter e and is 
in the same representation of the chiral symmetry as 2(27, i)- We perform global chiral and 
continuum fits using the form 

M^S=2 ^^^^ij^ 4'^°^°^ a^) + cithi + Chinih -mh^)+ c^m^ + Cy{my - ih]^) , (77) 

where m/ and are the sea and valence light-quark masses, ihh and thy the corresponding 
strange-quark masses and m/,g is the physical bare strange quark mass. The coefficients Ca 
depend on the action as indicated. By performing the global fits, c|P^^^ can be determined and 
the size of the lattice artefacts can be determined. Using all our data we find that the artefacts 
are 12% in the SU(2) chiral limit and 18% at the physical quark masses. If we restrict the data 
to pions with masses less than 350 Me V, we find artefacts of 10% in the chiral limit and 14% 
for physical quark masses. 

Based on these calculations we estimate the uncertainty due to the lattice artefacts as being 
15%, which we combine with the remaining uncertainties in quadrature. This estimate of the 
discretization error includes possible artefacts in the conversion of the renormalization constants 
from the IDSDR to the Iwasaki lattices. We stress that while lattice artefacts are the dominant 
source of systematic uncertainty in the present work, they will be reliably reduced when the 



45 





mi = 0.004 


mi = 0.006 


mi = 0.008 


Re(A2)xlO^ GeV 
Im(A2) X 10" GeV 


0.697(44) 
-14.73(37) 


0.748(41) 
-14.99(35) 


0.719(38) 
-15.23(34) 



TABLE IX: The amplitude A2, computed on the Iwasaki ensembles, after extrapolation to physical kaon 
and pion masses. The two pions in the final state are at rest (up to finite-volume effects) and energy is 
not conserved in these amplitudes (see text). 

calculations are repeated at a second lattice spacing. 



VI. ESTIMATING THE ERROR DUE TO PARTIAL QUENCHING 

The calculations described in this paper were designed to have almost physical kinematics, 
i.e. the kaon and pions have masses which are close to their physical values. This is achieved 
however, by the sea and valence quark masses being different; the sea-quark masses are mf^ = 
0.001 and mf"" = 0.045 and the valence masses are mY^^""^^ = 0.0001 and mf^''''^ = 0.049. 
Although we do not expect the dependence on the sea-quark to be very significant, in this 
section we report on some studies to check this. We start by describing an investigation of the 
sea-quark mass dependence performed with the 32^ Iwasaki lattices and in Subsec.VIB we 
report on the results obtained by reweighting mf^ from 0.001 to the valence value of 0.0001. 
Note that as the bare mass decreases from 0.001 to 0.0001, m/ -|-mres decreases by a relatively 
smaller ratio, from 0.0028 to 0.0019. 



A. Sea-quark mass dependence on the 32^ Iwasaki lattices 

K — )• TtTt correlation functions were also computed on the 32^ x 64, = 16 Iwasaki 
lattices (a^^ = 2.285(29) GeV) with three different light sea-quark masses m^'^^ = 
0.004,0.006,0.008 [|56l 1571 . For each of the sea-quark masses, the correlation functions were 
calculated using several valence masses: m^^i'^"^'^ = 0.002,0.004,0.006,0.008,0.025,0.03. Pe- 
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riodic boundary conditions were used, so the pions have zero momentum, resulting in a decay 
which does not conserve energy. For each of the three sea-quark masses, a chiral extrapolation 
was performed over the valence masses to determine the K — )■ kk amplitudes corresponding to 
physical kaon and pion masses (for the strange quark in the kaon this was an interpolation). The 
results are summarised in Tab.HXl 

From the table we see that any dependence on the light sea-quark mass is small, and gener- 
ally within the statistical uncertainties. As an estimate of the uncertainty we take the standard 
deviation of the results obtained with the different sea light-quark masses; 3.5% for Re(A2) and 
1 .7% for Im(A2). Although the kinematics are different from those for the physical decay on the 
IDSDR lattice, we take this to be an estimate of the error due to partial quenching. The range 
of sea-quark masses on the Iwasaki lattice and the long length of the extrapolation suggest that 
this may be a conservative estimate. We do not attempt to estimate the error due to the partial 
quenching of the strange quark, but note that the deviation from unitarity in the strange-quark 
mass is relatively small (m^^^ = 0.045 compared to mj^''^"'^'^ = 0.049) . 



B. Reweighting the light sea quarks 

The technique of reweighting allows us to change the sea-quark masses a posteriori, i.e. after 
the generation of the configurations [58], albeit at a loss of statistical precision. It is commonly 
used to correct for any difference between the simulated and physical strange-quark masses, see 
for example [fTTll . Here we reweight the light-quark mass in order to investigate the effects of 
its partial quenching. 

The reweighting is performed in 30 increments from the simulated mass mf'^ = 0.001 down 
to a value of mf^ = 0.0001 which corresponds to the valence light-quark mass and the results are 
shown in Fig. 1 1 The rightmost point in Fig.[TTJa) shows the result for Re A2 before reweight- 



ing, while the remaining points show the results after reweighting to the mass indicated on the 
X-axis, ending with m^*^^ = 0.0001 for the leftmost point. Similarly Fig.jTTj^b) shows the effects 
of reweighting on ImA2. The final results after reweighting are shown in Tab. [X] where they 
are compared with the results before reweighting. In this table, for illustration of the effects of 
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1.35 



101 (rw) 



0.001 (Sim) 



oM(. 



101 (rw) 



0.001 (Sim) 



(a) Reweighting Re A2 (b) Reweighting Im A2 

FIG. 11: Reweighting A2 from mf^ = 0.001 to mf ^ = 0.0001. 





m/ =0.001 


mi = 0.0001 (reweighted) 


ReA2 
ImA2 


1.381(38) X 10-*^ GeV 
-6.54(15) X 10-" GeV 


1.367(65) X 10-^ GeV 
-6.91(23) X 10-" GeV 



TABLE X: A2 before and after reweighting. The quoted errors correspond to the statistical fluctuations 
in the correlation functions only. The statistical uncertainties in the determination of the lattice spacing 
and non-perturbative renormalization have been omitted here. 

reweighting, we only quote the statistical error from the correlation functions themselves; we do 
not include the statistical errors from the determination of the lattice spacing or renormalization 
or any of the systematic errors. 

Examining the figures, it can be seen that, as expected, the statistical errors on ReA2 and 
ImA2 grow. Table |X] shows that the real part of A2 remains unchanged whereas the cen- 
tral value of the imaginary part decreases by 5.7% which is more than the 1.7% estimated 



in Sec VI A which we take to be our main estimate of the error due to partial quenching. An 
alternative approach would be to eliminate the systematic error due to partial quenching by us- 
ing the reweighted values for our final results. In doing this the systematic errors on ReA2 and 
ImA2 are unchanged at 18% and 19% respectively. Using the reweighted values, we would 
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ReA2 


ImA2 


lattice artefacts 


15% 


15% 


finite-volume corrections 


6.0% 


6.5% 


partial quenching 


3.5% 


1.7% 


renormalization 


1.8% 


5.6% 


unphysical kinematics 


0.4% 


0.8% 


derivative of the phase shift 


0.97% 


0.97% 


Wilson coefficients 


6.6% 


6.6% 


Total 


18% 


19% 



TABLE XI: Systematic error budget for Re A2 and ImA2. 

obtain the following results for the complex amplitude A2: 
ReA2 = 1.367(70)s,at(246)systl0-^GeV, ImAz = -6.91(51)s,at(131)systlO-i3GeV. (78) 



The results of Eq. (78 1 should be compared with Eq. (25 ), and it is clear that the differences due 



to reweighting are well within the total error. 



VII. ERROR BUDGET 



The sources of systematic error in the calculation of ReA2 and ImA2 include those from 
lattice artefacts, finite-volume effects, partial quenching, the uncertainty in the non-perturbative 
renormalization, the unphysical kinematics used in the calculation, the determination of the 
derivative of the phase shift and the Wilson coefficients. Although some of these uncertainties 
have been estimated in previous sections (NPR in Sec TV lattice artefacts in Sec jV] and partial 



quenching in Sec VI), here we summarise the conclusions of sections IV VI and briefly discuss 
the remaining sources of uncertainty before finaUy combining them all into a total systematic 
error. The results can be found in Tab. 
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A. Lattice Artefacts 



The estimate of the systematic error due to lattice artefacts is described in SecjV] and was 



estimated to be 15%. Comparing this with the other errors in Tab. XT we see that lattice artefacts 
are the dominant source of systematic error. They would be very significantly reduced by 
repeating the calculation at a second value of the lattice spacing. 

B. Finite- Volume Corrections 

In order to estimate the systematic error due to the finite volume of the lattice, we use SU(3) 
finite-volume chiral perturbation theory, in which the loop-integrals in Feynman diagrams are 
replaced by discrete sums over the allowed momenta. Expressions for the M = 3/2 K ^ nn 
matrix elements, ^(27,1) = (^^^ |2(27,i)l^^) ^i^d ^(8,8) = (^^^12(8,8)1^*^) known to 
next-to-leading order in SU (3) chiral perturbation theory. Since in chiral perturbation theory 
to leading order there is a single M = 3/2 operator constructed from the Goldstone boson 
fields which transforms as the (8,8) representation, the estimates derived below are the same 
for 2(8,8) ^iid 2(8.8)inix- There is also a single operator at lowest order which transforms as the 
(27,1) representation. We will be considering the leading order terms (labelled by "LO") and 
leading (one-loop) logarithmic terms (labelled by "log"). The LO expressions are well known 
and can be found in [|59l and [|60l . For ^^'27 1) ^1- ™ [l59l , (where we have added 

logarithmic terms from {mj^ — m^)i-\oop by hand as necessitated by Eq. (25) and corrected a 
factor of 1//^ in equation (A2)), and for -^^gl) we use Eq. (E3) in [[6OII . 

We denote the finite-volume corrections to the logarithmic terms in ^(27,1) and ^(8,8) by 
A^^'°^ and A^^^°|^ respectively. We estimate the relative size of these corrections, by using 
the pion and kaon masses in our lattice calculation finding, 

A^/°f,. A^/°|. 

^ - 0.0597 and ^ = 0.0649 (79) 
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if we normalize to the leading order expressions of the matrix elements, and 

^ -0.0352 and = 0.0438 (80) 



•^(27,1) +-^(27,1 



if we normalize to the leading order plus leading logarithmic expressions. More details can be 
found in llSTll . 

Evidently the leading logarithmic terms make significant corrections to the leading order 
terms. To have confidence that the chiral perturbation theory is converging we should check 
the size of the next-to-leading-order terms, but as these have unknown coefficients we are un- 
able to make a numerical estimate. We therefore make a conservative estimate by taking the 
larger relative finite- volume correction of Eq. f79] ) and conclude that the (27,1) operator carries 
a 6.0% finite-volume correction and that the (8,8) operator carries a 6.5% finite-volume correc- 
tion. Since ReA2 is dominated by the (27,1) operator and ImA2 is dominated by the (8,8)mix 
operator, these are the percentage errors due to finite- volume effects we assign to ReA2 and 
ImA2 respectively. 

C. Partial Quenching 

The effects of partial quenching have been discussed in detail in section fVl} Here we simply 
remind the reader that we neglect any systematic error due to partial quenching of the heavy- 
quark and attribute a 3.5% error to ReA2 and a 1.7% error to ImA2 due to the partial quenching 
in the light-quark sector of this calculation. 

D. Uncertainties due to the Renormalization 

We consider two sources of systematic error from the calculation of the renormalization 
constants. The first is designed to take into account lattice artefacts of higher order than G{a^) in 
the continuum extrapolation of the step-scaling function using the Iwasaki lattices, as described 



in section IV B[ and corresponds to the second error in equation (67 1. This systematic error is 
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ReA2 xlO^ GeV 


ImA2 X 10" GeV 


RI-SMOM(,^,^) 
RI-SMOM(7^,7^) 


1.381(46) Stat (15) (NPR-sy s) 
1.362(44)stat(03)(NPR-sys) 


-6.54(46)stat(33)(NPR-sys) 
-6.35(34)stat(42)(NPR-sys) 



TABLE XII: ReA2 and ImA2 calculated in the two different schemes. 



estimated in the same way that the statistical NPR error on A2 is calculated, i.e. Eq. (22) is 



used, but in this case dZ denotes the systematic errors on the Z- factors. The resulting error is 



displayed in Tab. XII and is labelled NPR-sys. We find this to be a 1.1% effect for ReA2 and a 
5.0% effect for ImA2 (see the second row of the table). 

The second source of systematic error in the renormalization constants is due to the trun- 
cation error in the perturbative matching to the MS scheme and to G{a^') scaling errors since 
we only have one lattice spacing and the Z-factors in the different schemes need not approach 
the continuum limit along the same scaling trajectory. Following conversion to the MS scheme, 
the four intermediate NPR schemes described in Sec IV A| should give equivalent answers. We 
estimate the resulting systematic error by considering the spread in results when A2 is calculated 
in the RI-SM0M(7^, y^^) scheme and in the RI-SMOM(^, ^) scheme. 

The results for A2 in the RI-SM0M(7^, 7^) and RI-SMOM(^,^) schemes are presented in 



Tab. XII We observe a 1.4% spread for Re A2 and a 2.5% spread for ImA2. Combining the two 



sources of error in quadrature, we find a 1.8% error for ReA2 and a 5.6% error for ImA2. 



E. Uncertainties due to the Unphysical Kinematics 

When choosing the parameters of the simulation, including the quark masses, the coupling 
constant and even the volume, we aim to obtain physical kaon and pion masses and E-,nz = niK- 
Once the simulation has been performed, we naturally find that this is not quite the case (see 
Tab.|l]) and we now attempt to estimate the systematic error that these non-physical kinematics 
contribute to our calculation. 

In addition to the results from the current simulation, we have a large collection ofK—^ nn 
amplitudes calculated on quenched lattices with a variety of light and strange-quark masses 
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and pion momenta. We use the observed dependence of the amplitudes with the param- 
eters to estimate our uncertainty due to the unphysical kinematics. On the quenched lat- 
tices we have a total of 60 values for the K nn amplitudes, obtained with all combi- 
nations of am, = 0.0023, 0.0047, 0.0071, am, = 0.046, 0.062, 0.078, 0.094, 0.110 and with 
n = 0, 1,2 and 3, where n is the number of spatial directions in which antiperiodic boundary 
conditions are imposed, n parametrizes the pion momenta as briefly explained in Sec II C[ 



The procedure for estimating the systematic error due to non-physical kinematics uses these 
quenched amplitudes, extrapolating the results in ami and interpolating them in am, and n, first 
to physical kinematics, and then to the kinematics simulated on the IDSDR lattices. This proce- 
dure is described in detail in [[57l . and is very similar to the extrapolation procedure described 



in section VI A when computing the error due to partial quenching. The difference here is that 
we can now interpolate to the correct pion-momenta. This is achieved by fitting the two-pion 
energy as a function of n, and interpolating to find ^P^y^ the value of n which corresponds to 
the desired two-pion energy. This in turn allows the decay amplitude to be interpolated and 
evaluated at n"^^^^. 

For the extrapolation to physical kinematics we find from the quenched lattices: 

ReA2 = 2.25 X lO^^GeV, ImA2 = -13.45 x 10"^^ GeV, (81) 

while the extrapolation to m^, mx and E-nn simulated in this article gives 

ReA2 = 2.26x lO^^GeV, ImA2 = -13.56 x lO^^^GeV. (82) 

We take the percentage differences between the two extrapolations as a measure of the system- 
atic error due to simulating at non-physical kinematics, and find 0.4% for ReA2 and 0.8% for 
ImA2. 



F. Uncertainty in tlie Derivative of tlie Pliase Shift 



The derivative of the s-wave phase shift dd/dk appearing in the Lellouch-Liischer factor was 
found by evaluating the derivative of the phenomenological curve at the momentum simulated 
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in our lattice calculation. This was discussed in section III and illustrated in Fig.[6j Alternatively 
we could have made a crude estimate of the derivative by taking the slope of the straight line 
between the phase shift at 17.63 MeV and 196.8 MeV. (c.f. the results of Tab|ll]). We estimate the 
systematic error to be 0.97%, which we find by calculating the percentage difference between 
the final results as obtained by the two different approaches. Since the derivative of the phase- 



shift only contributes a small fraction to the Lellouch-Liischer factor (see Tab. Ill I it is not 
surprising that the corresponding error is negligible. We note also that the derivative of the 
phase-shift can be calculated directly using the method proposed in [611 . 



G. Uncertainties in tlie evaluation of tlie Wilson coefficients 



The Wilson coefficients, which are calculated in perturbation theory and hence are not part 
of our lattice computations, are a necessary ingredient in the determination of the amplitude A2. 
The values presented in Tab. [V] were calculated at next-to-leading order (NLO) following the 
procedure outlined in [[34ll . In this section we estimate the systematic error due to the truncation 
of perturbation theory. To this end we calculate the Wilson coefficients to leading order (LO), 
following the procedure in [34] and measure the effect this has on the final results for ReA2 
and Im A2. The LO contribution to the Wilson coefficients is defined according to the following 
procedure: 

1. A value is chosen for the A parameter of four-flavor QCD. In ref. [[341 a range of values 
from 215 MeV to 435 MeV was used. In this paper we use the value of 328 MeV, which 
is close to the value corresponding to as{Mz) = 0. 1 184 [[39l . 

2. In setting the initial conditions for the Wilson coefficients at the scale of the W mass, 
corrections of 0{a) and 0{as) are only included when they depend on the top-quark 
mass. This also applies when calculating the coefficients zi at the scale of the charm mass 
(Eq.(VII.17) in[[34[[). 

3. In the QCD running to lower energies the one-loop expressions for the anomalous dimen- 
sion matrix and /3 -function are used. In the presence of electromagnetic interactions, the 



54 





LO 


NLO 


ReA2 
ImA2 


1.289(42) xlO"^ GeV 
-6.11(36)xlO-" GeV 


1.381(46)xlO-^ GeV 
-6.54(46) X 10-13 GeV 



TABLE XIII: Re(A2) and Im(A2) as calculated with LO Wilson coefficients and NLO Wilson coeffi- 
cients. The errors quoted here represent the total statistical uncertainty. 

LO anomalous dimension matrix also includes the term — ri^^ 

An ' 

4. At leading order the Wilson coefficients are continuous when crossing quark-mass thresh- 
olds. 



Tab. XIII shows how the decay amplitude varies when the LO Wilson coefficients are used 
instead of the NLO Wilson coefficients. The error in A2 due to the truncation in the perturbative 
calculation of the Wilson coefficients is very conservatively estimated by taking the difference 
between the NLO result and the LO result, and calculating this as a percentage of the LO result. 
We find the error to be 7. 1% for Re A2 and 8.1% for ImA2. 



Vin. SUMMARY AND CONCLUSIONS 

In ref. [[11 and the present paper we have presented the results of the first ab initio calculation 
of the complex K — )■ {mz)i=2 decay amplitude A2 and our results can be found in Eq. ([T]). It is 
very encouraging that our result for ReA2 agrees with the known experimental value and we 
are also able to determine ImA2 for the first time. The calculation was made possible by the 
major theoretical advances and technical progress which has been achieved over many years as 
described in the text above. Much of the important particle physics phenomenology, including 
the description of the weak interactions of the quarks in terms of matrix elements of specified 
four-quark operators multiplied by Wilson coefficients, has been understood since the 1970's, 
even before the methods of lattice QCD had been invented. However, it was only after major 
advances in lattice techniques that this calculation has become possible. The good control of 
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chiral symmetry provided by the 5-dimensional domain wall formulation, the ability to trans- 
late from lattice to continuum normalization of operators using non-perturbative methods and 
the finite-volume techniques capable of creating the proper, interacting kk final state with the 
correct energy are all essential ingredients in the calculation presented here. In addition, im- 
provements in computational algorithms, and teraflops-scale computing resources, enable us to 
perform the simulations with nearly physical u and d quark masses. 

The error on our result is dominated by lattice artefacts due to the fact that the calculation 
was performed at a single, rather course, lattice spacing a^^ ~ 1.4GeV. The most important 
extension of the calculation of A2 is therefore to repeat it at different values of /3, or at least 
at a second lattice spacing, so that the discretization errors can be essentially eliminated by 
extrapolating to zero lattice spacing. In addition, since the methods to compute A2 are now well 
in hand, more refined calculations using a larger lattice volume and physical light-quark masses 
(for the sea quarks as well as the valence ones) should be possible. These enhancements to 
the calculation reported here are well within reach of the next generation of high performance 
computers and should reduce the errors on the result for A2 by nearly an order of magnitude. 

Much more challenging but of even greater interest is the application of these methods to 
the calculation of the complex 7 = amplitude Aq. The calculation of both Aq and A2 from 
first principles will allow a direct comparison of e' / £ with the experimental result, giving new 
sensitivity to the search for physics beyond the Standard Model. The computational framework 
presented here will also support the calculation of Aq. However, serious obstacles must be 
overcome. Much larger Monte Carlo samples will be required to remove the large statistical 
fluctuations remaining after the contribution of the vacuum state has been removed. The device 
of applying anti-periodic boundary conditions to a single quark field used in this paper cannot 
be used in the case of the 1 = nn state. More sophisticated boundary conditions mixing 
quarks and anti-quarks and an isospin rotation, the so called G-parity boundary conditions, 
must be used instead for both the valence and the sea quarks. Exploratory studies (3 suggest 
that obtaining adequate Monte Carlo statistics will be practical with the next generation of high 
performance computers and efforts are presently underway to develop the necessary boundary 
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conditions. We anticipate that a complete calculation of CP violation in K nn decay within 
the Standard Model will be achieved before the fiftieth anniversary of its original discovery. 
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